#ifndef CGF_UTIL
#define CGF_UTIL

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/math/quaternion.hpp>

namespace cgmath {

	static const double PI = 3.141592653589793238462643383279502884;

	template<class T>
	T deg2Rad(T angle);

	template<class T>
	T rad2Deg(T angle);

	/*
	Calculates the signed distance between a plane defined by the point p0
	and the normal vector n to the point "point"
	http://softsurfer.com/Archive/algorithm_0104/algorithm_0104.htm
	*/
	template<class T>
	T signedPlaneDist(boost::numeric::ublas::vector<T> & p0,
					boost::numeric::ublas::vector<T> & n,
					boost::numeric::ublas::vector<T> & point);

	/* Creates A matrix for rotation around the x-axis.
	Real Time rendering p. 28
	*/
	template<class T>
	boost::numeric::ublas::matrix<T> rotXMat(T angle);

	/* Creates A matrix for rotation around the y-axis.
	Real Time rendering p. 28
	*/
	template<class T>
	boost::numeric::ublas::matrix<T> rotYMat(T angle);

	/* Creates A matrix for rotation around the z-axis.
	Real Time rendering p. 28
	*/
	template<class T>
	boost::numeric::ublas::matrix<T> rotZMat(T angle);

	/* Computes the dot product of two vectors
	Described in "Real-Time Rendering se" p.48 */
	template<class T>
	boost::numeric::ublas::vector<T> cross(boost::numeric::ublas::vector<T> & u,
		boost::numeric::ublas::vector<T> & v);

	/* Given three orthonormal vectors r,s,t and an angle delta
	it computes the matrix for rotating around v 
	Described in "Real-Time Rendering se" p.42 */
	template<class T>
	boost::numeric::ublas::matrix<T> rotMatrix(boost::numeric::ublas::vector<T> & r,
		boost::numeric::ublas::vector<T> & s,
		boost::numeric::ublas::vector<T> & t,
		T delta);

	/* Copy the elements of the array data to the vector target.
	The length of data must be at least the lenght of the vector*/
	template<class T>
	void arrayToVector(T * data,
		boost::numeric::ublas::vector<T> & target);

	/* Copy the elements of the array data to the matrix target.
	data = {1,2,3,4;
	5,6,7,8;
	9,10,11,12;
	13,14,15,16}*/
	template<class T>
	void arrayToMatrix(T * data,
		boost::numeric::ublas::matrix<T> & target);
	/* inverse operation */
	template<class T>
	void matrixToArray(T * data,
		boost::numeric::ublas::matrix<T> & target);

	/* Converts the quaternion to the target matrix.
	Described in "Real-Time Rendering se" p.48 and 
	http://www.flipcode.com/documents/matrfaq.html#Q54 */
	template<class T>
	void quaternionToMatrix(const boost::math::quaternion<T> & q, 
		boost::numeric::ublas::matrix<T> & targetMatrix);

	/* Convert the matrix to a quaternion   
	Described in "Real-Time Rendering se" p.49 and 
	http://www.flipcode.com/documents/matrfaq.html#Q55 */
	template<class T>
	void matrixToQuaternion(const boost::numeric::ublas::matrix<T> & matrix,
		boost::math::quaternion<T> & q);

	/* Extract the euler angels form the matrix.
	Described in "Real-Time Rendering se" p.40 */
	template<class T>
	void EulerAngles(const boost::numeric::ublas::matrix<T> & matrix,
		T & h,
		T & r,
		T & p);

	template<class T>
	void matrixFromAngles(boost::numeric::ublas::matrix<T> & matrix,
		const T & h,
		const T & r,
		const T & p);

	template<class T>
	void translationFromMatrix(const boost::numeric::ublas::matrix<T> & matrix,
		T& trX,
		T& trY,
		T& trZ);

	template<class T>
	void translationToMatrix(boost::numeric::ublas::matrix<T> & matrix,
		const T& trX,
		const T& trY,
		const T& trZ);

	template<class T>
	void printMatrix(const boost::numeric::ublas::matrix<T> & matrix);

	template<class T>
	/* Real-Time Rendering p.487 */
	T bernstein(unsigned int n, unsigned int i, T t);

	/* http://www.brpreiss.com/books/opus4/html/page467.html */
	unsigned int binom(unsigned int n, unsigned int m);
}

#endif /*CGF_UTIL*/
